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Abstract 



o 

A simple and effective method for imaging through ground-level atmospheric turbulence. 

1 Introduction 

The problem of reconstructing an image from a sequence that is altered by ground-level atmospheric 
^ turbulence is still largely unsolved, due to the time- and space-dependent nature of the distortion. 

In fact what is commonly referred to as turbulence in imaging can be modeled, at least in first 
approximation, as the combined effect of (i) a time-dependent deformation of the image domain 
and (ii) a blur with an anisoplanatic point spread function. This suggests an approach to solve the 
problem by first correcting for the geometric distortion and later applying a deblurring algorithm. 

Such approach is followed by other authors; see, for example, [HI El EH [18l EH]- In several 
cases, however, the geometric correction is achieved by simply computing the temporal mean of 
the data sequence. This introduces additional blur, whose effects, in the case of strong turbulence, 
are particularly destructive of image features. In this report we present a novel geometric idea, 
which is the one of correcting for the domain distortion by first computing the temporal mean 
of the deformations of the images with respect to the first one, and then creating the centroid 
image by warping the first one via this average deformation. At a second stage the data (i.e. the 
\Q \ image sequence) is registered onto the centroid and then averaged in time; this produces a blurri 

image that is made sharp by nonlocal Total Variation (NL-TV) deconvolution [SJE2]- In order to 
formulate all of this mathematically we need a notion of deformation; with this in mind, we briefly 
summarize the concept of optical flow and the standard techniques for its recovery. 



2 Optical flow 

Roughly speaking, optical flow between two images is "a vector field that deforms one image into 
the other". More precisely, given an image sequence 7(x,y,t), (x,y) E E [0,T], the optical flow 
between the images 7(x,y,to) and I(x,y,to + St), (x,y) E ^i, (for fixed to and St) is a vector field 
f (x, y) = (u(x : y),v(x, y)) such that: l(x + u(x, y),y + v(x, y), to + St) ~ /(x, y, to), (x, y) E J1. 

There is a very large and consolidated literature on the recovery of optical flow; see, for ex- 
ample, [21 13 E3 E51 E2 and references therein. A good part of it relies on the data correlation 
constraint, i.e. on the assumption that the image brightness of a moving region essentially remains 
constant in time. That is, if the location of a point moves within the image frame according 
to the map t \-> (x(t),y(t)), then it is the case that ^/(^(t), y(t), t) = 0. The chain rule yields: 
M^^fycft^ff where the velocity field (u, v) := (dx/dt, dy/dt) is interpreted as the optical 
flow at time t and location (x,y). If we rewrite such equation as 

^ + ^_l_^- (1) 
dx dy dt 
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Figure 1: (a) The centroid in Euclidean space, (b) P\C (thick) is obtained by averaging vectors P\P\ 
(thin), (c) The "correction" Q<5 (thick) is obtained by averaging vectors Qpi (thin). 

which must be solved in u and v, we note that the recovery of optical flow may be regarded as an 
inverse problem [HE]. In fact while equation (pQ) is the transport equation [4] for the function /, it is 
parameters (u, v) that must be recovered (the function / is known from the data). Like most inverse 
problems optical flow recovery is ill-posed in that at each location y) G ft the solution of (jTJ) is a 
line in the (u, i;)-plane. Moreover, in the texture free regions (characterized by dl/dx = dl/dy = 0) 
the problem is even more underdetermined. Therefore further hypotheses must be introduced. 

A fast and effective method was introduced by Horn & Schunk [ 9J, [TO] : given an image sequence / 
the optical flow at a fixed time t is a vector field (u, v) : R 2 that minimizes the energy 



E[u, v] = { («g + «| + |) 2 + «|| Vuf + a|| V,|| 2 } dx dy; 



(2) 



in other words we have added the i7 1 -seminorm of the optical flow (u, v) as regularizing term. The 
Euler-Lagrange equations of ([2]) are the partial differential equations: 

81/ dl dl dl\ A „ 81/ dl dl dl\ 

^~ + + ^7 ~ aAu = °- ^- + -z-v + — J - aAv 0. 

ox V ax ay at J ay V ox ay at J 

One sees immediately that in a texture-free region, in which the first terms on the left-hand side of 
the above equations are zero, the resulting optic flow is a vector field with harmonic components, 
whose boundary values are determined by the image regions that do have texture. Note also that 
the parameter a > essentially determines the smoothness of the solution (u,v). 

Other optical flow recovery methods have emerged: most notably, Black & Anandan [3l [16] 
have developed a robust scheme whose output is sharp and allows one to detect the boundaries of 
moving objects. In our application the motion field is produced by atmospheric turbulence and 
is thus relatively smooth by its own nature: whence Horn & Schunk's algorithm, which has the 
further advantage of being computationally efficient, will be our method of choice. 

3 The Centroid 

In Euclidean geometry the centroid of a set of points {Pi, . . . , Pjy} whose coordinates are known is 
given by C = ^ YliLi ^ see Figure [T](a). If only the positions of the points with respect to one of 

them, e.g. without loss of generality Pi, are known, i.e. if only vectors Pii^, i = 2, . . . N are known, 
then we can still compute the position of C with respect to P\ by calculating 



i+i 



(3) 

see Figure HJ^b). Also, it is immediate to verify that J^ii Cii = 0, and that Q(3 = ^2iLi QP\ 
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for any point Q; see Figure [T](c) . 

These elementary considerations lead to the following idea. Suppose we want to compute the 
"centroid" of the set of N images {Ii(x),x E &}i<i<N- Thinking of images as "points" and of 
the optical flow between two images as a "vector" between them, we can compute the optical 
flow between the reference image I\ and the rest of them, i.e. the N "vectors" I\fi(x), x E fi, 
i = 1, . . . , TV; so, treating optical flow as an element of a linear space, we compute the average 
flow I\& using equation This guiding principle may be summarized as follows: 



Instead of computing the temporal mean of the images, 
we compute the temporal mean of the deformations 
of the images with respect to the first one. 



Once this average deformation (optical flow) 1\C is known, and we warp I\ via the flow 1\C to 
finally get the centroid image C of the data set. Note however that since we are dealing with images 
(and not points in Euclidean space) this procedure is not exact. Most notably, the resulting C 
actually depends quite strongly on the choice of the initial reference image: if certain features are 
missing in I\ (due to the effects of optical turbulence), they will also be absent in the centroid, 
since the latter is obtained by warping I\. Also, it is the case that J2iLi C^i is usually not equal 
to zero. In fact, a "better" result is obtained by implementing the two following procedures. 



A. Iteration. The following iterative algorithm can be implemented: 

1. a first estimate of the centroid is computed with the procedure described above, and we shall 
call this Q; 

2. the correction flow Q$ := YliLi Q^i * s calculated (typically this is a non-zero vector field), 
see Figure [Jjc) ; a new estimate of the centroid image C is computed by warping Q via the 
optical flow QC; 

3. the previous step is repeated until convergence is achieved. 

In fact the correction will never become zero, but its L 2 norm will decrease at each step. Our 
experience with our data sets is that such norm typically stabilizes after about 4 or 5 iterations. 



B. Registration of the data onto the Centroid, and averaging. The image that is obtained 
by the previous procedure is better than the initial estimate of the centroid (i.e. the one obtained 
with only one iteration), in that the geometric correction is improved. However, as we noted above, 
since the centroid is obtained by warping one particular image of the data set, say Ji, it will lack 
the features that the optical turbulence has eliminated in such frame. A procedure to alleviate this 
problem is the following: 

1. register each image I{ onto the centroid C; this can be done by computing the flow tpi := cfi 
and defining if'ix) '•— h{x + ^(x)), x E 

2. average the resulting registered sequence, 

Under the assumption that all features of interest appear in the majority of the data frames, the 
above procedure allows one to recover them (albeit in blurred form) in the resulting image C R . 



3 



4 Deblurring 



So far we have used the data sequence ii, . . . , ijv to produce an image C R where the geometric 
distortion has been rectified, but that still is blurry. In other words, we have reduced our problem 
to a deblurring problem, for which we employ known deconvolution techniques. 

Given the fact that in most of our images there are repeated patterns, our method of choice is 
nonlocal Total Variation deconvolution (NL-TV); see O [12]. The energy functional is of the type 



E[u] := / (/ — k * u) 2 dx + a \ [u(x) — u(y)] 2 Wf(x,y) dy dx (4) 

where the weight is defined as 

, v / {Ga*\f(x + .)-f(y + .)\)(0 y 
Wf(x, y) := exp ^ — 



with x,y E ft. In the above formula G a is a Gaussian mask with standard deviation a and h 
is a scaling parameter. Note that the weight is approximately equal to 1 if the patches around 
locations x and y are similar, otherwise it is approximately zero. By the second term on the right- 
hand side of (j4j), that the weight "enforces" the same similarities that appear in the data / onto 
the minimizers u of the above energy. 

In our case the data / is the blurry image C R , whereas the k is a Guassian kernel whose standard 
deviation is chosen manually in a way that yields the most visually appealing result. In the future 
we plan to test physics-inspired point-spread functions, such as the Fried kernel described in [8j. 



5 Results 

In Figure E] we present the results for three data sets. For each one we show: 

(A) a frame from the original data sequence; 

(B) the temporal mean; 

(C) the centroid C; 

(D) the mean C R of the images obtained by registering the data set onto the centroid; 

(E) the image resulting from applying NL-TV to C R . 

In fact we never use the temporal mean in our algorithm, but we show it here to emphasize 
how the centroid is generally much sharper and preserves more features than the temporal mean; 
this is especially true for data sets with intense turbulence, such as the one shown in the second 
column of Figure EJ We should note that while the centroid contains artifacts that are inherited 
by the particular image that is used to create it, this in not the case with the image C R obtained 
by averaging the registered data. This is particularly apparent for the NATO sequence, shown in 
column (3) in Figure [2j note in particular how the second bar from the right has much more regular 
edges in image (D)-(3) than in image (C)-(3), and how some of the dots are better defined. 



6 Current and future work 

The geometric method that we have briefly described here proves quite effective for many data 
sets, such as the three that we have shown in this report. It is conceptually very simple and rather 
straightforward to implement. The results are comparable with the state of the art. 

Perhaps the most fundamental issue that is currently being explored is the dependence of the 
centroid on the initial (reference) image, as we already discussed above; the procedure described 
at the end of section [3] (denoted with B) certainly alleviates the problem as it allows some of the 
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possibly missing features to be recovered, but experimental results show that there is still a certain 
degree of dependence on the reference image. We are currently exploring algorithms to deal with 
this issue. Future work will also focus on speeding up the algorithm and selecting the parameters 
(e.g. in the computation of optical flow, or in the deconvolution step) in an automated manner. 
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